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^ . 

On . The sound radiation of 3 MHz acoustically driven air bubbles in 

liquid is analysed with respect to possible applications in second har- 
^1 monic ultrasound diagnostics devices, which have recently come into 

' clinical use. In the forcing pressure amplitude = 1 — 10 atm and 

■ ambient radius Rq = 0.5 — 5 fim parameter domain a narrow regime 
around the resonance radius Rq ^ 1 — 1.5 /im and relatively modest 
Pa ~ 2 — 2.5 atm is identified in which optimal sound yield in the second 

[ harmonic is achieved while maintaining spherical stability of the bubble. 

^T) I For smaller Pa and larger Rq hardly any sound is radiated; for larger 

■ Pa bubbles become unstable towards non-spherical shape oscillations of 

■ their surface. The computation of these instabilities is essential for the 
, evaluation of the optimal parameter regime. A region of slightly smaller 

Rq and Pq ~ 1 — 3 atm is best suited to achieve large ratios of the second 
harmonic to the fundamental intensity. Spherical stability is guaranteed 
' in the suggested regimes for liquids with an enhanced viscosity compared 



o 



to water, such as blood. 



I. INTRODUCTION 



^ I Microbubbles, i.e., gas bubbles of a few /im diameter, have long been known to 

d I be very effective scatterers of ultrasound (cf. e.g. Their scattering cross sec- 

tions for MHz sound waves can be more than two orders of magnitude greater than 
their geometrical cross sections 0. In the last decade, the concept of exploiting this 
property to perform refined ultrasound diagnostics with gas bubbles as echo contrast 
enhancers has enjoyed increasing attention . The general idea of this technique is to 
inject a microbubble suspension into a vein and to study the blood flow by detecting 
the bubbles' sound echo reaction to an applied acoustical field. This leads to ultra- 
sound images of higher contrast and quality as compared to conventional diagnostic 
techniques using only the ultrasound backscatter from the tissue itself. 

The quality of an ultrasonogram mainly depends on its spatial resolution and 
its signal intensity, more specifically, on the ratio of signal intensities from "desired" 
echoes (such as the blood flow in bubble diagnostics) to "background noise" reflections 
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(from surrounding tissue). Spatial resolution is, of course, limited by the wavelength 
of the ultrasound. But as the absorption of sound in tissue increases exponentially 
with frequency [^,|^ , frequencies below 10 MHz are used in most clinical applications. 
As a typical value, we will choose an ultrasound driving frequency of Udl^^n =3 MHz 
throughout this work. 

In doing diagnostics with bubble suspensions, it is desirable to improve the signal 
to noise ratio. Namely, when detecting the emitted sound from the bubble at the 
driving frequency 3 MHz, the signal is obscured by the driving and its reflections 
from tissue. To improve the signal quality, it has recently been proposed ||^,|^ to 
detect higher harmonics of the driving frequency in the sound emission spectrum of 
the bubble. In view of the aforementioned strong damping of higher frequencies, the 
lowest (second) harmonic at 6 MHz is of particular interest. It can be selectively 
excited if the parameters are chosen appropriately. We expect that soft tissue, driven 
into the regime of nonlinear response by the strong driving, will also reflect part of 
the sound in higher harmonics. Also, the large amplitude driving signal itself will 
undergo nonlinear distortion, generating higher harmonic frequency components [^. 
As in the case of the conventional method, it may therefore be necessary to focus on 
the difference between the reflected signal with and without injected microbubbles. 
However, by choosing the proper size of the bubbles and the proper forcing pressure 
amplitude, it is possible to enhance the bubbles' reflection signal in higher harmonics. 
The main focus of this paper is to suggest a parameter regime well suited for such an 
endeavor. 

Experimental and theoretical research on bubble dynamics has received consid- 
erable attention since the discovery of single bubble sonoluminescence by Gaitan in 
1990 (cf. 0) and the detailed experiments by the Putterman group at UCLA |[TD|-|T2|. 
In those experiments single microbubbles are driven with a frequency of ~ 30 kHz and 
with a pressure amplitude of Pa = 1.1 — 1.5 atm. Under very special conditions on 
experimental parameters such as the gas concentration in the liquid and the pressure 
amplitude, the emission of short light pulses (once per driving period) from the center 
of the bubble is observed. These experiments stimulated us to perform a series of 
studies on bubble stability ITSHT^. Three types of instability mechanisms seem to 
be important: Spherical instability, diffusive instability, and chemical instability. All 
of these studies are based on a Rayleigh-Plesset-like (RP) equation which provides 
an accurate description of the bubble wall dynamics even for strongly nonlinear os- 
cillations. Excellent agreement with the experiments was obtained, encouraging us 
to rely on the RP equation also for microbubbles driven at 3 MHz. We will give an 
overview on RP dynamics in section Section |ITl| shows results of calculated sound 
intensities emitted into the whole spectrum as well as at the fundamental and second 
harmonic frequencies. Of the instability mechanisms mentioned above, only shape 
instabilities are important here. They will be treated in section |^ and reveal im- 
portant restrictions on useful values of driving pressure amplitudes and bubble radii. 
Section presents conclusions. 
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II. RAYLEIGH-PLESSET BUBBLE DYNAMICS 



We will treat the dynamics and sound emission of a single bubble here; we assume 
that it is driven by a spatially homogeneous, standing wave field 

P{t) = PaCOSLUdt (1) 

with a frequency ujd/2n = 3MHz, corresponding to a period of T = 0.33 fis, and a 
sound amplitude Pa between 1 and lOatm (roughly 10^ to lO^Pa), refiecting typical 
peak pressures of devices in ultrasound diagnostics. Much higher pressure amplitudes, 



as applied in lithotripters (cf. e.g. |jT9|), could damage the tissue. 

Let us briefiy discuss the approximations we made here. In water, the chosen 
frequency corresponds to a sound wave length of A = 27iCi/uJd ~ 500 /im, where 
Q = 1481m/s is the sound velocity in water. The typical ambient radius Rq of the 
microbubbles is in the range of 1 — 5 fim. Therefore, the approximation of spatial 
homogeneity is justified. In diagnostic ultrasound devices, the driving sound is not a 
standing wave, but a short traveling wave pulse (which is not strictly monochromatic). 
The corresponding spatial fiuctuations of the driving pressure gradient at the location 
of the bubble will exert translational forces on the bubble. The resulting translational 
movements of the bubble are neglected (we will come back to this assumption later in 
this section) as well as bubble-bubble interactions, the so-called secondary Bjerknes 
forces [^. Finally, pressure fiuctuations due to the blood pressure (order of mag- 
nitude 0.05 atm) can also safely be neglected. In many cases, ultrasound contrast 
enhancers do not contain pure air bubbles, but stabilized bubbles with an albumin 
or saccharide coating This may lead to a shift in the resonance frequency of the 



bubbles as elaborated by de Jong, Church, and others [pl|-p^. 

Under these assumptions the dynamics of the bubble radius R{t) may be described 
by the following ordinary differential equation [^,^: 

RR +Ir^ = - {p{R, t) - P{t) - Po) 
^ Pi 

R d , , „ s „ / N N R 2cr , ^ 

+ — - {p{R, t) - Pit ) - Av- - — . 2 
Pici at R piR 

Typical parameters for an air bubble in water are the surface tension a = 0.073 kg/ s^, 
the water viscosity v = 10~^m^/s and density pi = 1000 kg /m^. We use these 
parameters for our calculations. Only the viscosity is chosen to be larger by a factor 
of three with respect to water (z/ = 3 • 10~^m^/s), corresponding to the value for 
blood. We will later see that the increased viscosity is essential for the spherical 
stability of the bubble at higher values of Pa- The external (ambient) pressure is 
Pq = 1 atm. We assume that the pressure inside the bubble is given by a van der 
Waals type equation of state 
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piRit)) = (Po + 



-)( 



Rl-h^ 



(3) 



with a (collective) van der Waals hard core radius h = i?o/8.54 (for air) i.e., 
is a measure for the total excluded volume of the molecules. The bubble radius under 
normal conditions (ambient radius) Rq is not uniform for the bubble population in 
a diagnostic suspension. The size distribution can, however, be controlled experi- 
mentally and is typically centered around 1-2 /im, with a width of about 1 /im [Q. 
For the effective polytropic exponent we take k ^ 1 as for the oscillation frequencies 
under consideration micrometer bubbles can be treated as approximately isothermal 
P7| . Equation (^ can be understood as a balance equation between the excitation 
due to the forcing (|l|) on the one hand and dissipative and acoustic loss processes on 
the other hand. 

In this paper, we will denote (^ the (modified) Rayleigh-Plesset (RP) equation, 
adopting a common practice in recent work on sonoluminescence [|^,^,^. Besides 
the pioneering work of Lord Rayleigh 
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and Plesset [^, other researchers have 
]) and a number of variations of this equation, e.g. Keller and Miksis 
or Gilmore 13311. Some of these variations are much more elaborate 



contributed to 

0, Flynn ^ _ 
than (^. However, a detailed comparison of the solutions obtained from these equa- 
tions |3^|l3 shows that significant deviations only occur for bubbles driven at very 
large pressure amplitudes ( > 5atm) and having large radii, i.e., Rq would have to be 
substantially larger than for the bubbles of the present study to necessitate the use 
of a more complicated dynamical equation. 
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FIG. 1. Potential according to eq. non-dimensionalized through dividing by ujd^, 
using Pa = 2 atm and Rq = 1.2 fim for three different phases w^t = 0, 7r/4, tt/2, respectively. 
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Let us first consider tlie resonance structure of tfie RP oscillator in the small 
forcing limit. To calculate the main resonance frequency we first note that for small 
forcing Pa < Pq the contribution of the sound coupling to the bubble dynamics (the 
term oc 1/q on the rhs of equation (^) is not important. In addition, if R stays large 
enough to ensure ^ (which is the case for weak driving), we can replace the 
van der Waals formula by an ideal gas expression. Writing R(t) = Ro{l + x{t)) we 
obtain 



X 



P, + -]{l + x)-'-^ Po + Pacosa;.t 2a 



RqJ 1 + X -Ro(l + x] 



Avx 3 X 



2 



Rlil + x) 21 + X 



(4) 



We interpret the bracketed term on the rhs of eq. (^) as an effective, time dependent 
force —dxV{x,t). Integration gives the time dependent potential 



V{x,t) 



^ fPo + ^) (1 + ^)""' + + Pa COSUdt) ln(l + x) - 

6K \ RqJ RqI + X 



(5) 



which displays a strong asymmetry in x, see Fig. |l]. If Pa = 0, the equilibrium point 
is a; = 0. For general Pa the minimum of the potential oscillates around this value. 
If Pa > Pq the potential is repulsive for a certain fraction of the period. For small 
Pa and thus small x we can linearize around x = and obtain a driven harmonic 
oscillator 

Pa 

X + 2'~^x + uJqX = — ^ cosudt (6) 

with the eigenfrequency 



^o = i/t1^(3/€Po + (3«:-1)— ) (7) 



2a \ 

piRl V^-"" ' '^-^ ^'RoJ 

and the damping constant 7 = 2v/R^. For fixed driving frequency ujd = 271 ■ SMHz 
and K = 1 the bubble oscillator is in (main) resonance if uod = uq. According to eq. 
(0), this corresponds to a resonance radius of Ro = 1.23 /zm. Taking viscous damping 
into account, the oscillation amplitude has its maximum value at a frequency 



^res 



which (with u = 3-10~^m^/s) shifts the main resonance radius to R^^^ = 1.18 fim. The 
corresponding radii for the subharmonics uj^H'^'' = uJd/2, uj^H'^^ = uJd/3, and toi^Jg'^^ = 
uJd/4: are Rq'^'^^ = 2.18 /xm, R'q^^^ = 3.14/im, and R'q^'^^ = 4.09 yum, respectively. The 
harmonic tu^g^ = 2uJd is at Rq'^ = 0.64 /xm. As we will see below, the harmonic and 
subharmonic resonances will also strongly affect the intensity of the emitted sound. 
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FIG. 2. The bubble radius R vs. time for four different parameter pairs {Ro,Pa); from 
upper to lower: a) {5 fim, 9.0 atm), b) (0.8 /xm, 1.2 atm), c) (1.2 ^um, 2.5 atm), d) (2.5 /im, 
4.5 atm) 

Figure 1^ shows the time series R{t) for four typical sets of parameters, as computed 
from (0) and (^ using a double precision, fourth-order, variable stepsize Runge-Kutta 
algorithm |^ . For small Pa it is trivial that the radius oscillates sinusoidally, but this 
can also happen for driving pressure amplitudes as large as = 9 atm, if the radius 
is much larger than the resonance radius, as seen in figure ^a). For other parameter 
combinations, the bubble changes its behavior and exhibits collapses, characterized 
by short duration minima of the radius, accompanied by large accelerations (large 
curvature of R(t)). In figure ||b) we observe one (weak) bubble collapse per cycle 
which becomes stronger for larger Pa, see figure ^c). For strong collapses the typical 
return time of the collapsing bubble is in the nanosecond range. Like the time series of 
most nonlinear oscillators, the bubble dynamics can period double so that a collapse 
only repeats every two cycles, as shown in figure |^d) for a strong collapse. In other 
parameter regions, aperiodic behavior ("chaos") can be observed (cf. |5^). It can 
also be seen from Fig. ^ that our notion of strong collapse coincides well with Flynn's 
32,38] definition of "transient cavities", for which the ratio of maximum expansion 
radius and ambient radius (expansion ratio) must fulfill Rmax/Ro ^ 2: the examples 
of Figs. I^c and d show rapid collapses and an expansion ratio of ~ 2. On the other 
hand. Figs. a and b exemplify weakly oscillating bubbles with expansion ratios near 



one. 



Rr 



For the purposes of this paper, it is instructive to compute the minimum radius 
lin = mmt{R(t)) which the bubble achieves during its oscillation. This quantity is 
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shown in Fig. |^. For this figure, as for the other 3D plots, the displayed function was 
evaluated at 100 x 100 equidistant grid points in the Pa — Rq plane. For weak forcing 
the minimum radius essentially equals the ambient radius (limit of small oscillations) . 
However, if the driving pressure amphtude is large enough, the bubble collapse is only 
halted in the immediate vicinity of the smallest possible radius, i.e., the van der Waals 
hard core radius. The transition towards hitting the hard core radius h = Rq/8.54: 
(upon increasing Pa) is rather abrupt, forming a well-defined threshold in the Pa - 
Rq plane. Obviously, this transition occurs for smaller Pa if the bubble radius is 
near one of the above mentioned resonance radii. The resonances at R^o^ = 1.18 /im 
and 4^/^^ = 2.18 IJ,m are clearly recognized in Fig. ^ As Pa becomes larger, the 



D 




FIG. 3. Minimum radius Rmin/Ro as a function of Rq and Pa- Note that the graph 
is enclosed between two planes: for small Pa, the quotient Rmin/Ro is nearly equal to one 
(weak oscillations), for large Pa and small Rq it approaches h/RQ = 1/8.54, with the van 
der Waals hard core radius h. Arrows in figures @, - ^ indicate axis orientation. 

We now come back to our above assumptions on the pressure field. The results 
of this work were obtained assuming a driving by standing plane waves. Today's 
ultrasound diagnostics devices usually emit a bundle of traveling waves that interfere 
constructively to build up large pressure peaks (5 — 10 atm) and to achieve sufficient 
spatial resolution. The resulting spatial pressure gradients lead to translational forces 
acting on the bubble, the so called primary Bjerknes forces [ pO|| : 
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FB{r,t) = V,{t)VP{r,t), 



(9) 



where Viy{t) is the time dependent volume of the bubble, and P{r,t) is the external 
pressure exerted on the bubble. However, numerical computation of shows that 
the accelerations resulting from a pressure gradient |VP| ~ Pa/^ are rather weak 
and that the translational velocities of the bubbles are small compared to their radial 
oscillation velocities. 



III. SOUND EMISSION OF OSCILLATORY BUBBLES 

Our focus of interest is on the sound emitted by the oscillating bubble. The 
far field sound pressure at a distance r ^ R from the center of the bubble can be 
calculated as [EIIBU 



An ultrasound diagnostics device will display a picture of sound intensity, which is 
based on the modulus (or, equivalently, the square) of the sound pressure. Obviously, 
the total detected intensity will not consist exclusively of signals due to (|10|) , but there 
will also be intensity components from reflections of the incoming signal (|l|) in the 
tissue. As these latter contributions depend on many peculiarities of the experimental 
or diagnostic setup, we do not try to model them here, but focus on the active sound 
radiation of the bubble. 

Figure ^ shows the total sound intensity / as a function of the forcing pressure 
amplitude Pa and the ambient radius Rq. According to Parseval's theorem it can 
be calculated either from the sound pressure time series Ps{r,t) or from its Fourier 
transform Ps{r,uj) = Ps{r,t) exp{iujt)dt, t ^T, namely 

1 rr 11 r+co 

I{r) = - \P,{r,t)\Ht = --— \P,{r,uo)\^duj. (11) 



T Jo T Zn J- 



oo 



We divide by the length r of the time series (r > 8T for all computations) and 
measure / in units of atm^. The intensity is evaluated at a distance of = 1cm 
from the bubble's center using standard double precision Fourier transform algorithms 
P^l . As / spans several orders of magnitude in our parameter regime, we also present 
its logarithm in Fig. Hb. 
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FIG. 4. a) Total sound intensity / at a distance of = 1cm from the bubble center 
as a function of Rq and Pa- The strong correlation between this figure and figure |3| is 
obvious. Note that the graph is truncated (higher /-values are not displayed) for better 
illustration of the resonance tongue structure. Figure b) shows a logarithmic plot of the 
same quantity. The small undulations at very large Pa on top of the resonance structure 
are due to numerical aliasing in the Fourier analysis and can be reduced with increasing 
computer power. 
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It can be seen from Fig. ^ that for small Pa or large Rq the bubble hardly emits 
any sound. Sound losses set in at a sharp threshold which is very similar to the 
threshold seen in Fig. |^ for the minimum radius. Moreover, comparing Figs. ^ and 
^ one realizes that strong sound emission and collapsing to the hard core radius 
are strongly correlated (note the opposite orientations of these two graphs). This 
behavior is expected from equation (|l3), as a stronger collapse means larger bubble 
wall acceleration R at the moment of collapse. The resonance structure in Rq is also 
clearly reflected in the emitted sound intensity. 
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FIG. 5. Time series of the sound pressure Ps{rN,t) from ( |To|) for the same four pairs of 
parameters as in figure ^. 

Time series of Ps{rN, t) and their power spectra for the same four parameter pairs 
{Rq, Pa) as in Fig. ^ are shown in Fig. ^ and Fig. ^ respectively. If the collapse 
is too violent, the sound intensity is distributed in a broad band spectrum with a 
large fraction of intensity emitted at higher harmonics, which will be absorbed by 
the tissue. Therefore it is appropriate to focus on the signal at the second harmonic 
frequency in order to get high intensities away from the driving frequency. The 
intensity hi'f') = \P{r, uj = 2uj(i)\'^ /t"^ in the second harmonic is displayed in Fig. |^. For 
comparison, we show the same plot also for the intensity Ii{r) = \P{r,uj = Udjl'^/T'^ 
of the fundamental (driving frequency) in Fig. ^. 
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FIG. 6. Power spectra \Ps{rN,i^d)\'^ for the same four pairs of parameters as in figure 




FIG. 7. Absolute intensity of the second harmonic of the sound emission (fQ). 
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Not surprisingly, the regions of greatest intensity of the fundamental as well as of 
the second harmonic coincide with those of the total intensity, i.e., they can be found 
at the resonance radii. The largest region and the highest maxima of second harmonic 
intensity occurs around R*^^ . For small Pa the intensity is of course nearly exclusively 
in the driving frequency itself, which is the frequency of the small oscillations of the 
bubbles. Upon increasing Pa, sound is also emitted in the second harmonic mode, 
for some parameter combinations up to 40% of the total intensity. At even larger 
Pa, higher and higher modes are excited, leaving a smaller and smaller fraction of 
total intensity for the second mode (cf. Fig. ^c and d). Even for large Pa, the 
driving frequency Ud remains the largest component of total emitted power in spite 
of the strong collapse, which displays much larger peak pressures, but only lasts for 
extremely short periods of time. 




FIG. 8. Absolute intensity of the fundamental of the sound emission (p^). 
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The key question now is: How should one choose the parameters for optimal 
detection of the second harmonic? The answer cannot be given exclusively from 
the absolute intensity of the second harmonic, Fig. |^. Clearly, the signal has to 
have a certain absolute intensity to overcome the noise level, and therefore Fig. |^ 
gives important information. One important application of the second harmonic 
method in ultrasound diagnostics, however, relies on the contrast between intensities 
at fundamental and second harmonic frequencies: When injecting a bubble suspension 
into the vascular system, a second harmonic signal is only expected from the contrast 
agent. Thus, detecting at 6 MHz in our example will give a bright image of the blood 
vessels. In a diagnostic situation, it should be possible to switch between this image 
and the scattering signal from surrounding tissue, which reflects the 3 MHz driving. 
But if the bubble emission signal at 3 MHz is more intense than these reflections, the 
vascular system will be the dominant feature in the fundamental frequency image, 
too. Therefore, meeting the demands of this application means to identify parameter 
regions where the ratio of second harmonic intensity to fundamental intensity /2//1 

is as high as possible. Fig. displays this quantity. It shows a distinct maximum at 

(2) 

small Rq ~ 0.7 — 0.8 fim close to the harmonic bubble resonance radius Rq . In this 
parameter region, the bubble essentially oscillates with 6 MHz instead of 3 MHz; an 
example for this behavior can be seen in Figs. Pd and Pd. For even smaller radii, a 
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rise in /2//1 indicates other resonances. All bubbles with high /2//1 ratios emit little 
absolute acoustic intensity (cf. Figs. |). 

We have presented a variety of intensity diagrams in Figs. ^ - ^ in order to meet 
the different demands of different experimental setups or diagnostic applications. 
Accordingly, the optimal parameter ranges for second harmonic sonography depend 
on the focus of interest: If the important quantity is relative intensity l^l^i-, one would 
pick the maximum of Fig. ^, i.e., bubbles with i?o ~ 0.5 — 1.2 /im and driving pressure 
amplitudes Pa ~ 1 — 3 atm. If absolute intensity is the key variable, one would choose 
the region around the main resonance radius in i?o, i.e., ~ 1.0 — 1.5 [im. Moreover, 
Fig. 1^ suggests choosing the pressure amplitude Pa as large as possible. However, as 
we will show in the next paragraph, bubble shape instabilities set an upper limit on 
practically useful Pa- Note that a change in the driving frequency would shift the 
resonance radii and, consequently, also the location of regions of maximum sound 
emission. If, for example, uj^ is smaller, the resonances are shifted towards larger i?o, 
see equation (|^). 

Also, the total amount of bubbles should be large enough to guarantee a strong 
signal, but low enough to prevent considerable bubble-bubble interaction. Therefore, 
it is of primary importance to assure that a high percentage of the generated bubbles 
is in the correct i?o regime. All other bubbles are essentially useless regarding the 
yield in the second harmonic and may even obscure the measurements. 



IV. SPHERICAL STABILITY 



To take advantage of the above suggested parameter regimes derived from RP dy- 
namics, the bubbles in these domains should be spherically stable, i.e., stable against 
the growth of non-spherical bubble deformations, which could eventually lead to bub- 
ble fragmentation and a breakdown of sound emission. The corresponding stability 
analysis has been performed in detail in |15[. We give a brief summary here. Consider 
a small distortion of the spherical interface R{t), 



i?(t)+a„(t)r„(0,, 



where Yn is a spherical (surface) harmonic of degree n. An approximate linearized 
dynamical equation of the distortion an{t) for each mode has been derived in 
following the pioneering work of Prosperetti |^ 



It reads (cf. || 



dn + Bn{t)dn - v4„(t)a„ = 



(12) 



with 

Anit) 



(n-l)- 

^ 'r p^R"" 



i?3 



(n - l)(n + 2) + 2n{n + 2)(n - 1) — 

R 



(13) 
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^ , , 3R 2u 



{n + 2){2n + 1) - 2n{n + 2)^- 



(14) 



Here, j3n = {n ~ l)(ri + l){n + 2) and 5 is a viscous boundary layer cutoff [ |15| |, 



,.,„,„(/Z«). (15) 

\\l Ud 2n) 

If the coefficients Anit) and -Bn(^) are periodic with period T, ([T2|) is an equation of 
Hill's type and instability occurs whenever the magnitude of the maximal eigenvalue 
of the Floquet transition matrix Fn{T) of eq. (|T2]) is larger than one. The Floquet 
transition matrix Fn{T) is defined by 



Fn{T) ^ . (16) 





Now for some parameter regimes the radius is not periodic with period T and thus 
the coefficients An{t) and Bn{t) are not either. Therefore, rather than calculating the 
Floquet matrix Fn{T) we calculate a transition matrix Fn{NT) with a large integer 
N (here N = 20) to average the dynamics. We numerically compute the eigenvalues 
of Fn{NT). The logarithm of the maximum eigenvalue can be understood as an 
approximate Lyapunov exponent. If it is positive, the mode a„(t) grows exponentially 
and the bubble is unstable towards the corresponding mode of shape oscillation. In 
Fig. 10 b and c we show the resulting stability diagrams for the second and the 



third mode (n = 2 and n = 3, respectively). Bubbles are shape unstable in the dark 
regions of the Pa — Rq plane, and shape stable in the white areas. Generally, the 
n = 2 mode is the most unstable one, but there are regimes where this does not hold. 

The most pronounced features of these stability diagrams are "tongues" of insta- 
bility. In the low Pa regime equation ([T2|) reduces to a Mathieu equation; in this case 
the tongues of instability are the well known Mathieu tongues, as demonstrated in 
pr3| . Large viscosity strongly damps out this tongue structure, as seen from com- 



paring the stability diagrams for water and blood (different viscosities) in Fig. p!0| a 
and |TU|b. In some regions of parameter space, stable and unstable points seem to 
be mixed erratically. This is due to long-periodic or chaotic bubble dynamics for 
these Pa — Rq combinations, for which the results of our stability analysis over 20T 
depend sensitively on the initial conditions, so that for slightly deviating parameters 
the stability behavior may be completely different. 
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FIG. 10. Stability diagram for the n = 2 mode for bubbles in water (a) and in blood 
(b) and for the n = 3 mode for bubbles in blood (c) according to the Floquet multipliers 
computed from equation (p^). In the white regions the bubbles are parametrically stable 
towards perturbations of the corresponding mode, in the dark regions they are parametri- 
cally unstable. Bubbles in water are much less stable than those in blood and higher modes 
Qn with n > 3 tend to be more stable than the second mode 02. We stress that details of 

as well as on the choice 



these stability diagrams depend on our approximation (12) - (|T^^ 
of the parameters of the liquid and the averaging time which is 20T here. Thus they should 
only be considered as a reflection of the general trend. 
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As we learn from this figure, bubbles in the Rq regime of maximum (absolute) yield 
in the second harmonic become spherically unstable for drivings with Pa > 2.5 atm. 
Below, the bubbles are stable in blood (due to its enhanced viscosity), whereas for 
water these bubbles are unstable with respect to the a2-mode even at Pa ~ 0.5 atm. 
Note that it is risky even to be close to an unstable regime, as the bubble may dif- 
fusionally shrink or grow into these regimes. This suggests that experiments with 
bubble suspensions in water will probably give misleading results (with respect to 
clinical applications) and should rather be carried out in blood or a fluid of cor- 
respondingly enhanced viscosity. In our calculations, we have assumed a viscosity 
of = 3 ■ lO^^m^/s, which is at the low end of typical measured blood viscosities 
(~ 3 — 5 • 10^^ m? / [0). Thus, we have determined lower bounds of the instability 
thresholds in Fig. and the actual regions of stability may be somewhat larger. 
Note also that stability will be enhanced when lower driving frequencies are used as 
this leads to a larger effective damping of surface oscillations (cf. ||TB|). 

Stability analysis shows that very large driving amplitudes (say 10 atm) will not 
help provide large response signals from the suspended bubbles. Instead, it may be 
useful to limit the amplitudes in the bubble regions to < 2.5 atm. 

The predictions for the region of optimal /2//1 intensity ratio remain virtually 
unaltered, as there is very little overlap of this region with the instability tongues. 
Indeed, according to the regime of large /2//1 identified above, shape instability in 
this regime is to be expected only in a tiny area of parameter space at Rq ~ 1.2 ^tti 
and Pa ~ 3 atm (see Fig. p!0|b). 

Besides spherical instability, diffusive instability and chemical instability also are 
matters of concern, as pointed out in detail for 26 kHz forced bubbles in refs. [pr5| , p!6| . 
Both types of instabilities will only be important on long timescales (milliseconds or 
longer) where our approximation of a stable standing wave is not appropriate anyhow. 
This is why we postpone the discussion of these instabilities to future work. 

V. SUMMARY AND CONCLUSIONS 

Bubble suspensions as contrast enhancers in ultrasound diagnostics are now state 
of the art. Improving the image quality by detecting the second harmonic of the 
driving frequency in the emitted sound spectrum is likely to lead to their further 
acceptance, as the advantages compared to conventional methods become even more 
pronounced. We have identified regions in parameter space, i.e., values for the driving 
pressure amplitude Pa and the ambient bubble radius Rq, where relatively high sound 
intensities at the second harmonic frequency are to be expected. These regions are 
intimately connected to the resonance structure of the bubble oscillator and to the 
collapse dynamics of the bubble. The best suited parameter regime to achieve a high 
absolute second harmonic intensity I2 is located around the main bubble resonance 
radius, i.e., Rq = 1.0 — 1.5 fim if a fixed driving frequency of 3 MHz is used. Requiring 
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bubble stability towards non-spherical perturbations limits useful driving pressures to 
a maximum of about 2.5 atm, if the bubbles oscillate in blood (cf. Fig. p!0|b). Here, we 
employed the previously specified values for p;, q, i/, cr to model the average properties 
of blood. Bubbles in this parameter range are stable in fluids with (at least) three 
times higher viscosity than water. 

The intensity ratio /2//1 is important for diagnostic purposes (switching between 
"background" and contrast agent images). It is optimized for bubbles around Rq ~ 
0.8 fim and Pa ~ 1 — 3 atm. Note that, again, there is an upper limit to the strength 
of optimal driving. 

We therefore suggest not to use very high (> 2.5 atm) pressure amplitudes; a 
gentler driving may lead to a better image quality. Also, the bubble radii should be 
somewhat smaller than those dominant in bubble suspensions used today (e.g. SH 
U 508 A with a radius distribution peak at ~ 1.4/im ||^), if maximum efficiency at 
3 MHz driving frequency is to be achieved. A narrower distribution around the peak 
(i.e., more uniform radii) would, of course, further amplify the sound signal. 
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